LOAD DATA

df <- readRDS("../data/models/social-risk-crash-rate-data.rds")

HOT-ENCODE YEAR

I will treat year as a factor and one-hot encode it.

df$year <- as.factor(df$year)
year_dummies <- model.matrix(~ year - 1, data = df)

df <- cbind(df[ , !(names(df) %in% c("year"))], year_dummies)

PREPARE TARGET

We will remove all possible target variables and keep only one per model training.

# Choose your target variable (e.g., crash rate per 1,000 residents)
target_var <- "crash_rate_per_1000"

# Remove all target variables except selected
cols_to_remove <- grep("per_1000", 
                       names(df), 
                       value = TRUE)
cols_to_remove <- setdiff(cols_to_remove, target_var) # keep this column

df <- df %>% select(-all_of(cols_to_remove),)

# Create feature matrix and target vector
X <- df %>% select(-target_var, -borough, -total_population, -geoid)
y <- df[[target_var]]

glimpse(X)
Rows: 13,518
Columns: 44
$ pct_male_population              <dbl> 12.460865, 11.694881, 12.229106, 1…
$ pct_white_population             <dbl> 4.1225202, 3.6815086, 2.5856754, 2…
$ pct_black_population             <dbl> 2.435429, 2.630197, 2.833673, 2.41…
$ pct_asian_population             <dbl> 0.19803330, 0.14388387, 0.23376782…
$ pct_hispanic_population          <dbl> 6.411328, 6.607147, 5.982424, 6.07…
$ pct_foreign_born                 <dbl> 2.909566, 2.442189, 2.187254, 2.20…
$ pct_age_under_18                 <dbl> 2.130762, 2.643626, 2.333613, 2.38…
$ pct_age_18_34                    <dbl> 1.715654, 1.653705, 1.882339, 1.96…
$ pct_age_35_64                    <dbl> 3.296112, 3.079115, 3.041014, 3.04…
$ pct_age_65_plus                  <dbl> 1.80895804, 1.78607845, 1.56929356…
$ median_income                    <dbl> 58582.658, 49964.513, 68000.000, 7…
$ pct_income_under_25k             <dbl> 0.5445916, 0.5544325, 0.5325841, 0…
$ pct_income_25k_75k               <dbl> 1.0568123, 1.1376418, 0.9533662, 0…
$ pct_below_poverty                <dbl> 1.9574830, 1.9510653, 1.8091597, 1…
$ median_gross_rent                <dbl> 1579.1133, 1524.3577, 1701.0000, 1…
$ pct_owner_occupied               <dbl> 1.33318303, 1.35699756, 1.58439699…
$ pct_renter_occupied              <dbl> 1.2085934, 1.2305140, 1.1518859, 1…
$ pct_less_than_hs                 <dbl> 1.4509748, 1.1951954, 1.1302166, 1…
$ pct_hs_diploma                   <dbl> 1.5576081, 1.4196542, 1.2704773, 1…
$ pct_some_college                 <dbl> 0.8473540, 0.8997538, 0.7256966, 0…
$ pct_associates_degree            <dbl> 0.4912749, 0.3280552, 0.3923234, 0…
$ pct_bachelors_degree             <dbl> 0.9482748, 0.9457966, 0.8537607, 0…
$ pct_graduate_degree              <dbl> 0.3998749, 0.7098271, 1.1037907, 0…
$ pct_in_labor_force               <dbl> 3.566504, 3.430191, 3.555304, 3.59…
$ pct_not_in_labor_force           <dbl> 3.448445, 3.326595, 3.162980, 3.10…
$ unemployment_rate                <dbl> 15.750133, 13.478747, 10.806175, 1…
$ pct_commute_short                <dbl> 0.7350082, 0.4604284, 0.1585556, 0…
$ pct_commute_medium               <dbl> 1.7061331, 1.6882374, 1.2521824, 1…
$ pct_commute_long                 <dbl> 2.477320, 2.685832, 3.553271, 3.73…
$ pct_carpool                      <dbl> 0.35798327, 0.31462606, 0.00000000…
$ pct_public_transit               <dbl> 2.056500, 2.540030, 2.898721, 2.36…
$ pct_walk                         <dbl> 0.37702494, 0.30695226, 0.13009688…
$ pct_bike                         <dbl> 0.00000000, 0.00000000, 0.00000000…
$ pct_work_from_home               <dbl> 0.01713750, 0.04604284, 0.04675356…
$ pct_vehicle                      <dbl> 2.0165122, 1.9222885, 2.1120415, 2…
$ post_pandemic                    <int> 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1…
$ poverty_vehicle_interaction      <dbl> 3.9472883, 3.7505104, 3.8210202, 3…
$ unemployment_vehicle_interaction <dbl> 31.760336, 25.910041, 22.823090, 2…
$ year2018                         <dbl> 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0…
$ year2019                         <dbl> 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0…
$ year2020                         <dbl> 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0…
$ year2021                         <dbl> 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0…
$ year2022                         <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0…
$ year2023                         <dbl> 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1…

HYPERPARAMETER TUNING

What This Does - Optuna (Python) handles the search space and Bayesian optimization. - The final best parameters are applied to fit the final_model. - Search space: Instead of predefined grids, trial$suggest_float() and trial$suggest_int() explore a range of values. - Best parameters: study$best_params holds the optimal hyperparameters.

## CONVERT TO DMATRIX
dtrain_all <- xgb.DMatrix(data = as.matrix(X), label = y)

## Start Python venv
reticulate::use_virtualenv("r-reticulate", required = TRUE)

## OPTUNA-BASED SPATIAL CV
optuna <- import("optuna")

boroughs <- unique(df$borough)
folds <- lapply(boroughs, function(b) which(df$borough != b))

# Optuna objective
objective <- function(trial) {
  params <- list(
    booster = "gbtree",
    eta = trial$suggest_float("eta", 0.01, 0.3, log = TRUE),
    max_depth = trial$suggest_int("max_depth", 3, 12),
    min_child_weight = trial$suggest_int("min_child_weight", 1, 10),
    subsample = trial$suggest_float("subsample", 0.5, 1.0),
    colsample_bytree = trial$suggest_float("colsample_bytree", 0.5, 1.0),
    gamma = trial$suggest_float("gamma", 0, 10),
    lambda = trial$suggest_float("lambda", 0, 10),
    alpha = trial$suggest_float("alpha", 0, 10)
  )
  
  rmse_scores <- numeric(length(folds))
  for (i in seq_along(folds)) {
    train_idx <- folds[[i]]
    valid_idx <- setdiff(seq_len(nrow(dtrain_all)), train_idx)

    dtrain <- xgb.DMatrix(data = as.matrix(X[train_idx, ]), label = y[train_idx])
    dvalid <- xgb.DMatrix(data = as.matrix(X[valid_idx, ]), label = y[valid_idx])

    model <- xgb.train(
      params = params,
      data = dtrain,
      nrounds = 500,
      watchlist = list(val = dvalid),
      early_stopping_rounds = 20,
      verbose = 0
    )

    rmse_scores[i] <- min(model$evaluation_log$val_rmse)
  }
  
  preds <- predict(model, as.matrix(X[valid_idx, ]))
  return(Metrics::rmse(y[valid_idx], preds))
}

# Run Optuna study
set.seed(2025)
study <- optuna$create_study(direction = "minimize")
study$optimize(objective, n_trials = 50)

best_params <- study$best_params
print(best_params)

$eta [1] 0.2309768

$max_depth [1] 4

$min_child_weight [1] 3

$subsample [1] 0.5583975

$colsample_bytree [1] 0.8797697

$gamma [1] 3.35839

$lambda [1] 6.131609

$alpha [1] 2.805485

TRAIN/TEST SPLIT

# Set seed
set.seed(2025)

# Split by index
train_index <- createDataPartition(y, p = 0.8, list = FALSE)
X_train <- X[train_index, ]
y_train <- y[train_index]
X_test <- X[-train_index, ]
y_test <- y[-train_index]

# Convert to xgb.DMatrix
dtrain <- xgb.DMatrix(data = as.matrix(X_train), label = y_train)
dtest <- xgb.DMatrix(data = as.matrix(X_test), label = y_test)

TRAIN MODEL

# Set seed
set.seed(2025)

# Training with parallel processing
final_model <- xgb.train(
  params = list(
    eta = best_params$eta,
    max_depth = best_params$max_depth,
    gamma = best_params$gamma,
    colsample_bytree = best_params$colsample_bytree,
    min_child_weight = best_params$min_child_weight,
    subsample = best_params$subsample,
    objective = "reg:squarederror",
    eval_metric = "rmse"
  ),
  data = dtrain,
  nrounds = 1000,
  watchlist = list(train = dtrain, test = dtest),
  early_stopping_rounds = 20,
  verbose = 1,
  nthread = detectCores() - 1
)

[78] train-rmse:1.220267 test-rmse:1.555561 [79] train-rmse:1.211549 test-rmse:1.556347 [80] train-rmse:1.207152 test-rmse:1.553078 [81] train-rmse:1.205943 test-rmse:1.553145 [82] train-rmse:1.205311 test-rmse:1.552420 [83] train-rmse:1.201741 test-rmse:1.555835 [84] train-rmse:1.197931 test-rmse:1.553410 [85] train-rmse:1.191889 test-rmse:1.551760 [86] train-rmse:1.189975 test-rmse:1.551783 [87] train-rmse:1.188499 test-rmse:1.550743 Stopping. Best iteration: [67] train-rmse:1.258358 test-rmse:1.549545 ## SAVE MODEL

# Create directory if it doesn't exist
if (!dir.exists("../data/models")) {
  dir.create("../data/models", recursive = TRUE)
}

# Save the final XGBoost model
saveRDS(final_model, file = "../data/models/xgb_model.rds")

# Save the best parameters
saveRDS(best_params, file = "../data/models/xgb_best_params.rds")

cat("Model and parameters saved to ../data/models/")
# Load the final XGBoost model
final_model <- readRDS("../data/models/xgb_model.rds")

# Load the best parameters
best_params <- readRDS("../data/models/xgb_best_params.rds")

MODEL EVALUATION

library(Metrics)
library(ggplot2)
library(dplyr)

set.seed(2025)

# Predict on test set
preds <- predict(final_model, as.matrix(X_test))

# --- Metrics ---
rmse <- sqrt(mean((y_test - preds)^2))
mae <- mean(abs(y_test - preds))
mape <- mean(abs((y_test - preds) / y_test)) * 100
r2 <- 1 - (sum((y_test - preds)^2) / sum((y_test - mean(y_test))^2))

cat("Model Evaluation Metrics:\n")
Model Evaluation Metrics:
cat("  RMSE:", rmse, "\n")
  RMSE: 1.549545 
cat("  MAE :", mae, "\n")
  MAE : 0.8372956 
cat("  MAPE:", mape, "%\n")
  MAPE: Inf %
cat("  R²  :", r2, "\n\n")
  R²  : 0.2595503 
# --- Residuals ---
residuals <- y_test - preds
residual_df <- data.frame(
  actual = y_test,
  predicted = preds,
  residuals = residuals
)

# --- Plot: Predicted vs Actual ---
p1 <- residual_df %>%
  ggplot(aes(x = actual, y = predicted)) +
  geom_point(alpha = 0.5) +
  geom_abline(slope = 1, intercept = 0, color = "red") +
  theme_minimal() +
  labs(title = "Predicted vs Actual Crash Rates",
       x = "Actual",
       y = "Predicted")

# --- Plot: Residuals vs Predicted ---
p2 <- residual_df %>%
  ggplot(aes(x = predicted, y = residuals)) +
  geom_point(alpha = 0.5, color = "blue") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  theme_minimal() +
  labs(title = "Residuals vs Predicted",
       x = "Predicted",
       y = "Residuals")

# --- Plot: Residual Density ---
# Residual Histogram
p3 <- ggplot(residual_df, aes(x = residuals)) +
    geom_histogram(binwidth = 0.2, fill = "steelblue", color = "white") +
    geom_density(color = "red") +
    theme_minimal() +
    labs(title = "Residual Distribution", x = "Residuals", y = "Count")

# Print plots
print(p1)

print(p2)

print(p3)

ggsave("../report/plots/predicted_vs_actual_values_plot.png", p1, width = 10, height = 6, dpi = 300)
ggsave("../report/plots/resisuals_vs_predicted_values_plot.png", p2, width = 10, height = 6, dpi = 300)
ggsave("../report/plots/residual_density_plot.png", p3, width = 10, height = 6, dpi = 300)

Model Evaluation Metrics: RMSE: 1.549545 MAE : 0.8372956 MAPE: Inf % R² : 0.2595503 ## SHAP EXPLAINABILITY

# Compute SHAP values
shap_values <- shap.values(xgb_model = final_model, X_train = as.matrix(X_train))
shap_long <- shap.prep(shap_contrib = shap_values$shap_score, X_train = as.matrix(X_train))

# SHAP summary plot
print(shap.plot.summary(shap_long))

if (!dir.exists("../report/plots")) {
  dir.create("../report/plots")
}

shap <- shap.plot.summary(shap_long)
ggsave("../report/plots/shap_summary_plot.png", shap, width = 10, height = 6, dpi = 300)

INSPECT TREE ENSEMBLE

xgb.plot.tree(model = final_model, trees = 0)
xgb.plot.tree(model = final_model, trees = 1)
xgb.plot.tree(model = final_model, trees = 2)
xgb.plot.multi.trees(model = final_model)
# ============================================================
# Additional Model Diagnostics and Deeper Analysis
# ============================================================

library(ggplot2)
library(dplyr)
library(pdp)        # For Partial Dependence Plots
library(DALEX)      # For model explainability
library(ggthemes)
library(sf)

# ---------------------------
# 1. SHAP Dependence and Interaction Plots
# ---------------------------
message("\nGenerating SHAP dependence and interaction plots...")

# Assuming shap_values and shap_long are already computed
# (if not, recompute them using iml or SHAPforxgboost packages)

# Top feature by SHAP importance
top_feature <- shap_long %>%
  as_tibble() %>%
  count(variable, wt = abs(value), sort = TRUE) %>%
  dplyr::slice(1) %>%
  pull(variable)

# Dependence plot for top feature
shap.plot.dependence(data_long = shap_long, x = top_feature, color_feature = top_feature)


# Interaction values
shap_interaction_values <- predict(
  final_model,
  as.matrix(X_train),
  predinteraction = TRUE
)

# shap_interaction_values will be a 3D array: [n_samples, n_features, n_features]
interactions <- dim(shap_interaction_values)
# ---------------------------
# 3. Partial Dependence Plots (PDP)
# ---------------------------
message("\nGenerating Partial Dependence Plots...")

top_features <- shap_long %>%
  count(variable, wt = abs(value), sort = TRUE) %>%
  dplyr::slice(1:10) %>%
  pull(variable)

# Ensure the output directory exists
dir.create("../report/plots/pdp", recursive = TRUE, showWarnings = FALSE)

for (f in top_features) {
  pd <- partial(final_model, pred.var = f, train = as.matrix(X_train), grid.resolution = 30)
  p <- plot(pd, main = paste("Partial Dependence of", f))
  
  # Save as PNG
  ggsave(
    filename = paste0("../report/plots/pdp/pdp_", f, ".png"),
    plot = p,
    width = 8,
    height = 6,
    dpi = 300
  )
}

LS0tCnRpdGxlOiAiWEdCb29zdCBNb2RlbGluZyBSIE5vdGVib29rIgphdXRob3I6ICJBSiBTdHJhdW1hbi1TY290dCIKZGF0ZTogImByIFN5cy5EYXRlKClgIgpvdXRwdXQ6CiAgaHRtbF9kb2N1bWVudDoKICAgIGRmX3ByaW50OiBwYWdlZAogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQKICBwZGZfZG9jdW1lbnQ6IGRlZmF1bHQKLS0tCgpgYGB7ciBzZXR1cCwgaW5jbHVkZT1GQUxTRX0Ka25pdHI6Om9wdHNfY2h1bmskc2V0KGVjaG8gPSBUUlVFLCB3YXJuaW5nID0gRkFMU0UsIG1lc3NhZ2UgPSBGQUxTRSkKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkoY2FyZXQpCmxpYnJhcnkoeGdib29zdCkKbGlicmFyeShTSEFQZm9yeGdib29zdCkKbGlicmFyeShNYXRyaXgpCmxpYnJhcnkocmV0aWN1bGF0ZSkKbGlicmFyeShkb1BhcmFsbGVsKQpgYGAKCiMjIExPQUQgREFUQQoKYGBge3IgbG9hZC1kYXRhfQpkZiA8LSByZWFkUkRTKCIuLi9kYXRhL21vZGVscy9zb2NpYWwtcmlzay1jcmFzaC1yYXRlLWRhdGEucmRzIikKYGBgCgojIyBIT1QtRU5DT0RFIGBZRUFSYAoKSSB3aWxsIHRyZWF0IGB5ZWFyYCBhcyBhIGZhY3RvciBhbmQgb25lLWhvdCBlbmNvZGUgaXQuCgpgYGB7ciBob3QtZW5jb2RlLXllYXJ9CmRmJHllYXIgPC0gYXMuZmFjdG9yKGRmJHllYXIpCnllYXJfZHVtbWllcyA8LSBtb2RlbC5tYXRyaXgofiB5ZWFyIC0gMSwgZGF0YSA9IGRmKQoKZGYgPC0gY2JpbmQoZGZbICwgIShuYW1lcyhkZikgJWluJSBjKCJ5ZWFyIikpXSwgeWVhcl9kdW1taWVzKQpgYGAKCiMjIFBSRVBBUkUgVEFSR0VUCgpXZSB3aWxsIHJlbW92ZSBhbGwgcG9zc2libGUgdGFyZ2V0IHZhcmlhYmxlcyBhbmQga2VlcCBvbmx5IG9uZSBwZXIgbW9kZWwgdHJhaW5pbmcuCgpgYGB7ciBzaW5nbGUtdGFyZ2V0LWRmfQojIENob29zZSB5b3VyIHRhcmdldCB2YXJpYWJsZSAoZS5nLiwgY3Jhc2ggcmF0ZSBwZXIgMSwwMDAgcmVzaWRlbnRzKQp0YXJnZXRfdmFyIDwtICJjcmFzaF9yYXRlX3Blcl8xMDAwIgoKIyBSZW1vdmUgYWxsIHRhcmdldCB2YXJpYWJsZXMgZXhjZXB0IHNlbGVjdGVkCmNvbHNfdG9fcmVtb3ZlIDwtIGdyZXAoInBlcl8xMDAwIiwgCiAgICAgICAgICAgICAgICAgICAgICAgbmFtZXMoZGYpLCAKICAgICAgICAgICAgICAgICAgICAgICB2YWx1ZSA9IFRSVUUpCmNvbHNfdG9fcmVtb3ZlIDwtIHNldGRpZmYoY29sc190b19yZW1vdmUsIHRhcmdldF92YXIpICMga2VlcCB0aGlzIGNvbHVtbgoKZGYgPC0gZGYgJT4lIHNlbGVjdCgtYWxsX29mKGNvbHNfdG9fcmVtb3ZlKSwpCgojIENyZWF0ZSBmZWF0dXJlIG1hdHJpeCBhbmQgdGFyZ2V0IHZlY3RvcgpYIDwtIGRmICU+JSBzZWxlY3QoLXRhcmdldF92YXIsIC1ib3JvdWdoLCAtdG90YWxfcG9wdWxhdGlvbiwgLWdlb2lkKQp5IDwtIGRmW1t0YXJnZXRfdmFyXV0KCmdsaW1wc2UoWCkKYGBgCgoKIyMgSFlQRVJQQVJBTUVURVIgVFVOSU5HCgoqKldoYXQgVGhpcyBEb2VzKioKLSBPcHR1bmEgKFB5dGhvbikgaGFuZGxlcyB0aGUgc2VhcmNoIHNwYWNlIGFuZCBCYXllc2lhbiBvcHRpbWl6YXRpb24uCi0gVGhlIGZpbmFsIGJlc3QgcGFyYW1ldGVycyBhcmUgYXBwbGllZCB0byBmaXQgdGhlIGBmaW5hbF9tb2RlbGAuCi0gKipTZWFyY2ggc3BhY2U6KiogSW5zdGVhZCBvZiBwcmVkZWZpbmVkIGdyaWRzLCBgdHJpYWwkc3VnZ2VzdF9mbG9hdCgpYCBhbmQgYHRyaWFsJHN1Z2dlc3RfaW50KClgIGV4cGxvcmUgYSByYW5nZSBvZiB2YWx1ZXMuCi0gKipCZXN0IHBhcmFtZXRlcnM6KiogYHN0dWR5JGJlc3RfcGFyYW1zYCBob2xkcyB0aGUgb3B0aW1hbCBoeXBlcnBhcmFtZXRlcnMuCgpgYGB7ciB0dW5lLXBhcmFtcywgZXZhbD1GQUxTRX0KIyMgQ09OVkVSVCBUTyBETUFUUklYCmR0cmFpbl9hbGwgPC0geGdiLkRNYXRyaXgoZGF0YSA9IGFzLm1hdHJpeChYKSwgbGFiZWwgPSB5KQoKIyMgU3RhcnQgUHl0aG9uIHZlbnYKcmV0aWN1bGF0ZTo6dXNlX3ZpcnR1YWxlbnYoInItcmV0aWN1bGF0ZSIsIHJlcXVpcmVkID0gVFJVRSkKCiMjIE9QVFVOQS1CQVNFRCBTUEFUSUFMIENWCm9wdHVuYSA8LSBpbXBvcnQoIm9wdHVuYSIpCgpib3JvdWdocyA8LSB1bmlxdWUoZGYkYm9yb3VnaCkKZm9sZHMgPC0gbGFwcGx5KGJvcm91Z2hzLCBmdW5jdGlvbihiKSB3aGljaChkZiRib3JvdWdoICE9IGIpKQoKIyBPcHR1bmEgb2JqZWN0aXZlCm9iamVjdGl2ZSA8LSBmdW5jdGlvbih0cmlhbCkgewogIHBhcmFtcyA8LSBsaXN0KAogICAgYm9vc3RlciA9ICJnYnRyZWUiLAogICAgZXRhID0gdHJpYWwkc3VnZ2VzdF9mbG9hdCgiZXRhIiwgMC4wMSwgMC4zLCBsb2cgPSBUUlVFKSwKICAgIG1heF9kZXB0aCA9IHRyaWFsJHN1Z2dlc3RfaW50KCJtYXhfZGVwdGgiLCAzLCAxMiksCiAgICBtaW5fY2hpbGRfd2VpZ2h0ID0gdHJpYWwkc3VnZ2VzdF9pbnQoIm1pbl9jaGlsZF93ZWlnaHQiLCAxLCAxMCksCiAgICBzdWJzYW1wbGUgPSB0cmlhbCRzdWdnZXN0X2Zsb2F0KCJzdWJzYW1wbGUiLCAwLjUsIDEuMCksCiAgICBjb2xzYW1wbGVfYnl0cmVlID0gdHJpYWwkc3VnZ2VzdF9mbG9hdCgiY29sc2FtcGxlX2J5dHJlZSIsIDAuNSwgMS4wKSwKICAgIGdhbW1hID0gdHJpYWwkc3VnZ2VzdF9mbG9hdCgiZ2FtbWEiLCAwLCAxMCksCiAgICBsYW1iZGEgPSB0cmlhbCRzdWdnZXN0X2Zsb2F0KCJsYW1iZGEiLCAwLCAxMCksCiAgICBhbHBoYSA9IHRyaWFsJHN1Z2dlc3RfZmxvYXQoImFscGhhIiwgMCwgMTApCiAgKQogIAogIHJtc2Vfc2NvcmVzIDwtIG51bWVyaWMobGVuZ3RoKGZvbGRzKSkKICBmb3IgKGkgaW4gc2VxX2Fsb25nKGZvbGRzKSkgewogICAgdHJhaW5faWR4IDwtIGZvbGRzW1tpXV0KICAgIHZhbGlkX2lkeCA8LSBzZXRkaWZmKHNlcV9sZW4obnJvdyhkdHJhaW5fYWxsKSksIHRyYWluX2lkeCkKCiAgICBkdHJhaW4gPC0geGdiLkRNYXRyaXgoZGF0YSA9IGFzLm1hdHJpeChYW3RyYWluX2lkeCwgXSksIGxhYmVsID0geVt0cmFpbl9pZHhdKQogICAgZHZhbGlkIDwtIHhnYi5ETWF0cml4KGRhdGEgPSBhcy5tYXRyaXgoWFt2YWxpZF9pZHgsIF0pLCBsYWJlbCA9IHlbdmFsaWRfaWR4XSkKCiAgICBtb2RlbCA8LSB4Z2IudHJhaW4oCiAgICAgIHBhcmFtcyA9IHBhcmFtcywKICAgICAgZGF0YSA9IGR0cmFpbiwKICAgICAgbnJvdW5kcyA9IDUwMCwKICAgICAgd2F0Y2hsaXN0ID0gbGlzdCh2YWwgPSBkdmFsaWQpLAogICAgICBlYXJseV9zdG9wcGluZ19yb3VuZHMgPSAyMCwKICAgICAgdmVyYm9zZSA9IDAKICAgICkKCiAgICBybXNlX3Njb3Jlc1tpXSA8LSBtaW4obW9kZWwkZXZhbHVhdGlvbl9sb2ckdmFsX3Jtc2UpCiAgfQogIAogIHByZWRzIDwtIHByZWRpY3QobW9kZWwsIGFzLm1hdHJpeChYW3ZhbGlkX2lkeCwgXSkpCiAgcmV0dXJuKE1ldHJpY3M6OnJtc2UoeVt2YWxpZF9pZHhdLCBwcmVkcykpCn0KCiMgUnVuIE9wdHVuYSBzdHVkeQpzZXQuc2VlZCgyMDI1KQpzdHVkeSA8LSBvcHR1bmEkY3JlYXRlX3N0dWR5KGRpcmVjdGlvbiA9ICJtaW5pbWl6ZSIpCnN0dWR5JG9wdGltaXplKG9iamVjdGl2ZSwgbl90cmlhbHMgPSA1MCkKCmJlc3RfcGFyYW1zIDwtIHN0dWR5JGJlc3RfcGFyYW1zCnByaW50KGJlc3RfcGFyYW1zKQpgYGAKJGV0YQpbMV0gMC4yMzA5NzY4CgokbWF4X2RlcHRoClsxXSA0CgokbWluX2NoaWxkX3dlaWdodApbMV0gMwoKJHN1YnNhbXBsZQpbMV0gMC41NTgzOTc1CgokY29sc2FtcGxlX2J5dHJlZQpbMV0gMC44Nzk3Njk3CgokZ2FtbWEKWzFdIDMuMzU4MzkKCiRsYW1iZGEKWzFdIDYuMTMxNjA5CgokYWxwaGEKWzFdIDIuODA1NDg1CgojIyBUUkFJTi9URVNUIFNQTElUCgpgYGB7ciB0cmFpbi10ZXN0LXNwbGl0fQojIFNldCBzZWVkCnNldC5zZWVkKDIwMjUpCgojIFNwbGl0IGJ5IGluZGV4CnRyYWluX2luZGV4IDwtIGNyZWF0ZURhdGFQYXJ0aXRpb24oeSwgcCA9IDAuOCwgbGlzdCA9IEZBTFNFKQpYX3RyYWluIDwtIFhbdHJhaW5faW5kZXgsIF0KeV90cmFpbiA8LSB5W3RyYWluX2luZGV4XQpYX3Rlc3QgPC0gWFstdHJhaW5faW5kZXgsIF0KeV90ZXN0IDwtIHlbLXRyYWluX2luZGV4XQoKIyBDb252ZXJ0IHRvIHhnYi5ETWF0cml4CmR0cmFpbiA8LSB4Z2IuRE1hdHJpeChkYXRhID0gYXMubWF0cml4KFhfdHJhaW4pLCBsYWJlbCA9IHlfdHJhaW4pCmR0ZXN0IDwtIHhnYi5ETWF0cml4KGRhdGEgPSBhcy5tYXRyaXgoWF90ZXN0KSwgbGFiZWwgPSB5X3Rlc3QpCmBgYAoKIyMgVFJBSU4gTU9ERUwKCmBgYHtyIHRyYWluLW1vZGVsLCBldmFsPUZBTFNFfQojIFNldCBzZWVkCnNldC5zZWVkKDIwMjUpCgojIFRyYWluaW5nIHdpdGggcGFyYWxsZWwgcHJvY2Vzc2luZwpmaW5hbF9tb2RlbCA8LSB4Z2IudHJhaW4oCiAgcGFyYW1zID0gbGlzdCgKICAgIGV0YSA9IGJlc3RfcGFyYW1zJGV0YSwKICAgIG1heF9kZXB0aCA9IGJlc3RfcGFyYW1zJG1heF9kZXB0aCwKICAgIGdhbW1hID0gYmVzdF9wYXJhbXMkZ2FtbWEsCiAgICBjb2xzYW1wbGVfYnl0cmVlID0gYmVzdF9wYXJhbXMkY29sc2FtcGxlX2J5dHJlZSwKICAgIG1pbl9jaGlsZF93ZWlnaHQgPSBiZXN0X3BhcmFtcyRtaW5fY2hpbGRfd2VpZ2h0LAogICAgc3Vic2FtcGxlID0gYmVzdF9wYXJhbXMkc3Vic2FtcGxlLAogICAgb2JqZWN0aXZlID0gInJlZzpzcXVhcmVkZXJyb3IiLAogICAgZXZhbF9tZXRyaWMgPSAicm1zZSIKICApLAogIGRhdGEgPSBkdHJhaW4sCiAgbnJvdW5kcyA9IDEwMDAsCiAgd2F0Y2hsaXN0ID0gbGlzdCh0cmFpbiA9IGR0cmFpbiwgdGVzdCA9IGR0ZXN0KSwKICBlYXJseV9zdG9wcGluZ19yb3VuZHMgPSAyMCwKICB2ZXJib3NlID0gMSwKICBudGhyZWFkID0gZGV0ZWN0Q29yZXMoKSAtIDEKKQpgYGAKWzc4XQl0cmFpbi1ybXNlOjEuMjIwMjY3CXRlc3Qtcm1zZToxLjU1NTU2MSAKWzc5XQl0cmFpbi1ybXNlOjEuMjExNTQ5CXRlc3Qtcm1zZToxLjU1NjM0NyAKWzgwXQl0cmFpbi1ybXNlOjEuMjA3MTUyCXRlc3Qtcm1zZToxLjU1MzA3OCAKWzgxXQl0cmFpbi1ybXNlOjEuMjA1OTQzCXRlc3Qtcm1zZToxLjU1MzE0NSAKWzgyXQl0cmFpbi1ybXNlOjEuMjA1MzExCXRlc3Qtcm1zZToxLjU1MjQyMCAKWzgzXQl0cmFpbi1ybXNlOjEuMjAxNzQxCXRlc3Qtcm1zZToxLjU1NTgzNSAKWzg0XQl0cmFpbi1ybXNlOjEuMTk3OTMxCXRlc3Qtcm1zZToxLjU1MzQxMCAKWzg1XQl0cmFpbi1ybXNlOjEuMTkxODg5CXRlc3Qtcm1zZToxLjU1MTc2MCAKWzg2XQl0cmFpbi1ybXNlOjEuMTg5OTc1CXRlc3Qtcm1zZToxLjU1MTc4MyAKWzg3XQl0cmFpbi1ybXNlOjEuMTg4NDk5CXRlc3Qtcm1zZToxLjU1MDc0MyAKU3RvcHBpbmcuIEJlc3QgaXRlcmF0aW9uOgpbNjddCXRyYWluLXJtc2U6MS4yNTgzNTgJdGVzdC1ybXNlOjEuNTQ5NTQ1CiMjIFNBVkUgTU9ERUwKCmBgYHtyIHNhdmUtbW9kZWwsIGV2YWw9RkFMU0V9CiMgQ3JlYXRlIGRpcmVjdG9yeSBpZiBpdCBkb2Vzbid0IGV4aXN0CmlmICghZGlyLmV4aXN0cygiLi4vZGF0YS9tb2RlbHMiKSkgewogIGRpci5jcmVhdGUoIi4uL2RhdGEvbW9kZWxzIiwgcmVjdXJzaXZlID0gVFJVRSkKfQoKIyBTYXZlIHRoZSBmaW5hbCBYR0Jvb3N0IG1vZGVsCnNhdmVSRFMoZmluYWxfbW9kZWwsIGZpbGUgPSAiLi4vZGF0YS9tb2RlbHMveGdiX21vZGVsLnJkcyIpCgojIFNhdmUgdGhlIGJlc3QgcGFyYW1ldGVycwpzYXZlUkRTKGJlc3RfcGFyYW1zLCBmaWxlID0gIi4uL2RhdGEvbW9kZWxzL3hnYl9iZXN0X3BhcmFtcy5yZHMiKQoKY2F0KCJNb2RlbCBhbmQgcGFyYW1ldGVycyBzYXZlZCB0byAuLi9kYXRhL21vZGVscy8iKQpgYGAKCmBgYHtyIGxvYWQtbW9kZWx9CiMgTG9hZCB0aGUgZmluYWwgWEdCb29zdCBtb2RlbApmaW5hbF9tb2RlbCA8LSByZWFkUkRTKCIuLi9kYXRhL21vZGVscy94Z2JfbW9kZWwucmRzIikKCiMgTG9hZCB0aGUgYmVzdCBwYXJhbWV0ZXJzCmJlc3RfcGFyYW1zIDwtIHJlYWRSRFMoIi4uL2RhdGEvbW9kZWxzL3hnYl9iZXN0X3BhcmFtcy5yZHMiKQoKYGBgCgojIyBNT0RFTCBFVkFMVUFUSU9OCgpgYGB7ciBldmFsLW1vZGVsLCBtZXNzYWdlPUZBTFNFLCB3YXJuaW5nPUZBTFNFfQpsaWJyYXJ5KE1ldHJpY3MpCmxpYnJhcnkoZ2dwbG90MikKbGlicmFyeShkcGx5cikKCnNldC5zZWVkKDIwMjUpCgojIFByZWRpY3Qgb24gdGVzdCBzZXQKcHJlZHMgPC0gcHJlZGljdChmaW5hbF9tb2RlbCwgYXMubWF0cml4KFhfdGVzdCkpCgojIC0tLSBNZXRyaWNzIC0tLQpybXNlIDwtIHNxcnQobWVhbigoeV90ZXN0IC0gcHJlZHMpXjIpKQptYWUgPC0gbWVhbihhYnMoeV90ZXN0IC0gcHJlZHMpKQptYXBlIDwtIG1lYW4oYWJzKCh5X3Rlc3QgLSBwcmVkcykgLyB5X3Rlc3QpKSAqIDEwMApyMiA8LSAxIC0gKHN1bSgoeV90ZXN0IC0gcHJlZHMpXjIpIC8gc3VtKCh5X3Rlc3QgLSBtZWFuKHlfdGVzdCkpXjIpKQoKY2F0KCJNb2RlbCBFdmFsdWF0aW9uIE1ldHJpY3M6XG4iKQpjYXQoIiAgUk1TRToiLCBybXNlLCAiXG4iKQpjYXQoIiAgTUFFIDoiLCBtYWUsICJcbiIpCmNhdCgiICBNQVBFOiIsIG1hcGUsICIlXG4iKQpjYXQoIiAgUsKyICA6IiwgcjIsICJcblxuIikKCiMgLS0tIFJlc2lkdWFscyAtLS0KcmVzaWR1YWxzIDwtIHlfdGVzdCAtIHByZWRzCnJlc2lkdWFsX2RmIDwtIGRhdGEuZnJhbWUoCiAgYWN0dWFsID0geV90ZXN0LAogIHByZWRpY3RlZCA9IHByZWRzLAogIHJlc2lkdWFscyA9IHJlc2lkdWFscwopCgojIC0tLSBQbG90OiBQcmVkaWN0ZWQgdnMgQWN0dWFsIC0tLQpwMSA8LSByZXNpZHVhbF9kZiAlPiUKICBnZ3Bsb3QoYWVzKHggPSBhY3R1YWwsIHkgPSBwcmVkaWN0ZWQpKSArCiAgZ2VvbV9wb2ludChhbHBoYSA9IDAuNSkgKwogIGdlb21fYWJsaW5lKHNsb3BlID0gMSwgaW50ZXJjZXB0ID0gMCwgY29sb3IgPSAicmVkIikgKwogIHRoZW1lX21pbmltYWwoKSArCiAgbGFicyh0aXRsZSA9ICJQcmVkaWN0ZWQgdnMgQWN0dWFsIENyYXNoIFJhdGVzIiwKICAgICAgIHggPSAiQWN0dWFsIiwKICAgICAgIHkgPSAiUHJlZGljdGVkIikKCiMgLS0tIFBsb3Q6IFJlc2lkdWFscyB2cyBQcmVkaWN0ZWQgLS0tCnAyIDwtIHJlc2lkdWFsX2RmICU+JQogIGdncGxvdChhZXMoeCA9IHByZWRpY3RlZCwgeSA9IHJlc2lkdWFscykpICsKICBnZW9tX3BvaW50KGFscGhhID0gMC41LCBjb2xvciA9ICJibHVlIikgKwogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IDAsIGxpbmV0eXBlID0gImRhc2hlZCIsIGNvbG9yID0gInJlZCIpICsKICB0aGVtZV9taW5pbWFsKCkgKwogIGxhYnModGl0bGUgPSAiUmVzaWR1YWxzIHZzIFByZWRpY3RlZCIsCiAgICAgICB4ID0gIlByZWRpY3RlZCIsCiAgICAgICB5ID0gIlJlc2lkdWFscyIpCgojIC0tLSBQbG90OiBSZXNpZHVhbCBEZW5zaXR5IC0tLQojIFJlc2lkdWFsIEhpc3RvZ3JhbQpwMyA8LSBnZ3Bsb3QocmVzaWR1YWxfZGYsIGFlcyh4ID0gcmVzaWR1YWxzKSkgKwogICAgZ2VvbV9oaXN0b2dyYW0oYmlud2lkdGggPSAwLjIsIGZpbGwgPSAic3RlZWxibHVlIiwgY29sb3IgPSAid2hpdGUiKSArCiAgICBnZW9tX2RlbnNpdHkoY29sb3IgPSAicmVkIikgKwogICAgdGhlbWVfbWluaW1hbCgpICsKICAgIGxhYnModGl0bGUgPSAiUmVzaWR1YWwgRGlzdHJpYnV0aW9uIiwgeCA9ICJSZXNpZHVhbHMiLCB5ID0gIkNvdW50IikKCiMgUHJpbnQgcGxvdHMKcHJpbnQocDEpCnByaW50KHAyKQpwcmludChwMykKCmdnc2F2ZSgiLi4vcmVwb3J0L3Bsb3RzL3ByZWRpY3RlZF92c19hY3R1YWxfdmFsdWVzX3Bsb3QucG5nIiwgcDEsIHdpZHRoID0gMTAsIGhlaWdodCA9IDYsIGRwaSA9IDMwMCkKZ2dzYXZlKCIuLi9yZXBvcnQvcGxvdHMvcmVzaXN1YWxzX3ZzX3ByZWRpY3RlZF92YWx1ZXNfcGxvdC5wbmciLCBwMiwgd2lkdGggPSAxMCwgaGVpZ2h0ID0gNiwgZHBpID0gMzAwKQpnZ3NhdmUoIi4uL3JlcG9ydC9wbG90cy9yZXNpZHVhbF9kZW5zaXR5X3Bsb3QucG5nIiwgcDMsIHdpZHRoID0gMTAsIGhlaWdodCA9IDYsIGRwaSA9IDMwMCkKYGBgCk1vZGVsIEV2YWx1YXRpb24gTWV0cmljczoKICBSTVNFOiAxLjU0OTU0NSAKICBNQUUgOiAwLjgzNzI5NTYgCiAgTUFQRTogSW5mICUKICBSwrIgIDogMC4yNTk1NTAzIAojIyBTSEFQIEVYUExBSU5BQklMSVRZCgpgYGB7ciBzaGFwfQojIENvbXB1dGUgU0hBUCB2YWx1ZXMKc2hhcF92YWx1ZXMgPC0gc2hhcC52YWx1ZXMoeGdiX21vZGVsID0gZmluYWxfbW9kZWwsIFhfdHJhaW4gPSBhcy5tYXRyaXgoWF90cmFpbikpCnNoYXBfbG9uZyA8LSBzaGFwLnByZXAoc2hhcF9jb250cmliID0gc2hhcF92YWx1ZXMkc2hhcF9zY29yZSwgWF90cmFpbiA9IGFzLm1hdHJpeChYX3RyYWluKSkKCiMgU0hBUCBzdW1tYXJ5IHBsb3QKcHJpbnQoc2hhcC5wbG90LnN1bW1hcnkoc2hhcF9sb25nKSkKCmlmICghZGlyLmV4aXN0cygiLi4vcmVwb3J0L3Bsb3RzIikpIHsKICBkaXIuY3JlYXRlKCIuLi9yZXBvcnQvcGxvdHMiKQp9CgpzaGFwIDwtIHNoYXAucGxvdC5zdW1tYXJ5KHNoYXBfbG9uZykKZ2dzYXZlKCIuLi9yZXBvcnQvcGxvdHMvc2hhcF9zdW1tYXJ5X3Bsb3QucG5nIiwgc2hhcCwgd2lkdGggPSAxMCwgaGVpZ2h0ID0gNiwgZHBpID0gMzAwKQpgYGAKCiMjIElOU1BFQ1QgVFJFRSBFTlNFTUJMRQoKYGBge3IgaW5zcGVjdC10cmVlfQp4Z2IucGxvdC50cmVlKG1vZGVsID0gZmluYWxfbW9kZWwsIHRyZWVzID0gMCkKYGBgCmBgYHtyfQp4Z2IucGxvdC50cmVlKG1vZGVsID0gZmluYWxfbW9kZWwsIHRyZWVzID0gMSkKYGBgCmBgYHtyfQp4Z2IucGxvdC50cmVlKG1vZGVsID0gZmluYWxfbW9kZWwsIHRyZWVzID0gMikKYGBgCmBgYHtyfQp4Z2IucGxvdC5tdWx0aS50cmVlcyhtb2RlbCA9IGZpbmFsX21vZGVsKQpgYGAKCmBgYHtyfQojID09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PQojIEFkZGl0aW9uYWwgTW9kZWwgRGlhZ25vc3RpY3MgYW5kIERlZXBlciBBbmFseXNpcwojID09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PQoKbGlicmFyeShnZ3Bsb3QyKQpsaWJyYXJ5KGRwbHlyKQpsaWJyYXJ5KHBkcCkgICAgICAgICMgRm9yIFBhcnRpYWwgRGVwZW5kZW5jZSBQbG90cwpsaWJyYXJ5KERBTEVYKSAgICAgICMgRm9yIG1vZGVsIGV4cGxhaW5hYmlsaXR5CmxpYnJhcnkoZ2d0aGVtZXMpCmxpYnJhcnkoc2YpCgojIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQojIDEuIFNIQVAgRGVwZW5kZW5jZSBhbmQgSW50ZXJhY3Rpb24gUGxvdHMKIyAtLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0KbWVzc2FnZSgiXG5HZW5lcmF0aW5nIFNIQVAgZGVwZW5kZW5jZSBhbmQgaW50ZXJhY3Rpb24gcGxvdHMuLi4iKQoKIyBBc3N1bWluZyBzaGFwX3ZhbHVlcyBhbmQgc2hhcF9sb25nIGFyZSBhbHJlYWR5IGNvbXB1dGVkCiMgKGlmIG5vdCwgcmVjb21wdXRlIHRoZW0gdXNpbmcgaW1sIG9yIFNIQVBmb3J4Z2Jvb3N0IHBhY2thZ2VzKQoKIyBUb3AgZmVhdHVyZSBieSBTSEFQIGltcG9ydGFuY2UKdG9wX2ZlYXR1cmUgPC0gc2hhcF9sb25nICU+JQogIGFzX3RpYmJsZSgpICU+JQogIGNvdW50KHZhcmlhYmxlLCB3dCA9IGFicyh2YWx1ZSksIHNvcnQgPSBUUlVFKSAlPiUKICBkcGx5cjo6c2xpY2UoMSkgJT4lCiAgcHVsbCh2YXJpYWJsZSkKCiMgRGVwZW5kZW5jZSBwbG90IGZvciB0b3AgZmVhdHVyZQpzaGFwLnBsb3QuZGVwZW5kZW5jZShkYXRhX2xvbmcgPSBzaGFwX2xvbmcsIHggPSB0b3BfZmVhdHVyZSwgY29sb3JfZmVhdHVyZSA9IHRvcF9mZWF0dXJlKQoKIyBJbnRlcmFjdGlvbiB2YWx1ZXMKc2hhcF9pbnRlcmFjdGlvbl92YWx1ZXMgPC0gcHJlZGljdCgKICBmaW5hbF9tb2RlbCwKICBhcy5tYXRyaXgoWF90cmFpbiksCiAgcHJlZGludGVyYWN0aW9uID0gVFJVRQopCgojIHNoYXBfaW50ZXJhY3Rpb25fdmFsdWVzIHdpbGwgYmUgYSAzRCBhcnJheTogW25fc2FtcGxlcywgbl9mZWF0dXJlcywgbl9mZWF0dXJlc10KaW50ZXJhY3Rpb25zIDwtIGRpbShzaGFwX2ludGVyYWN0aW9uX3ZhbHVlcykKc2F2ZVJEUyhpbnRlcmFjdGlvbnMsICIuLi9kYXRhL21vZGVscy9zaGFwX2ludGVyYWN0aW9ucy5yZHMiKQpgYGAKCmBgYHtyfQojIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQojIDMuIFBhcnRpYWwgRGVwZW5kZW5jZSBQbG90cyAoUERQKQojIC0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLS0tLQptZXNzYWdlKCJcbkdlbmVyYXRpbmcgUGFydGlhbCBEZXBlbmRlbmNlIFBsb3RzLi4uIikKCnRvcF9mZWF0dXJlcyA8LSBzaGFwX2xvbmcgJT4lCiAgY291bnQodmFyaWFibGUsIHd0ID0gYWJzKHZhbHVlKSwgc29ydCA9IFRSVUUpICU+JQogIGRwbHlyOjpzbGljZSgxOjEwKSAlPiUKICBwdWxsKHZhcmlhYmxlKQoKIyBFbnN1cmUgdGhlIG91dHB1dCBkaXJlY3RvcnkgZXhpc3RzCmRpci5jcmVhdGUoIi4uL3JlcG9ydC9wbG90cy9wZHAiLCByZWN1cnNpdmUgPSBUUlVFLCBzaG93V2FybmluZ3MgPSBGQUxTRSkKCmZvciAoZiBpbiB0b3BfZmVhdHVyZXMpIHsKICBwZCA8LSBwYXJ0aWFsKGZpbmFsX21vZGVsLCBwcmVkLnZhciA9IGYsIHRyYWluID0gYXMubWF0cml4KFhfdHJhaW4pLCBncmlkLnJlc29sdXRpb24gPSAzMCkKICBwIDwtIHBsb3QocGQsIG1haW4gPSBwYXN0ZSgiUGFydGlhbCBEZXBlbmRlbmNlIG9mIiwgZikpCiAgCiAgIyBTYXZlIGFzIFBORwogIGdnc2F2ZSgKICAgIGZpbGVuYW1lID0gcGFzdGUwKCIuLi9yZXBvcnQvcGxvdHMvcGRwL3BkcF8iLCBmLCAiLnBuZyIpLAogICAgcGxvdCA9IHAsCiAgICB3aWR0aCA9IDgsCiAgICBoZWlnaHQgPSA2LAogICAgZHBpID0gMzAwCiAgKQp9CmBgYAo=